C MODPATH Version 3.00 (V3, Release 2, 5-99)
C Changes:
C   No change from previous release: (V3, Release 1, 9-94)
C***** SUBROUTINES *****
C     FLOLIN
C     NEWXYZ
C     DTCALC
C     ISIZS
C***********************
 
C***** SUBROUTINE *****
      SUBROUTINE FLOLIN (IPART,IGRID,MODE,TIME,TMAX,
     1IDSCH,JP,IP,KP,XP,YP,ZP,ZLOC,IBOUND,LAYCON,ZBOT,ZTOP,XMAX,YMAX,
     2QX,QY,QZ,DELC,DELR,POR,HEAD,NCON,NCOL,NROW,NLAY,NLPOR,NZDIM,
     3NCP1,NRP1,NLP1,ISNK,IREV,FRAC,IZSTOP,I2,I7,QSS,QSTO,INILOC,
     4NPART,ISS,ICASE,HDRY,NSTEP,ICMPCT)
C
      DIMENSION IBOUND(NCOL,NROW,NLAY),ZTOP(NZDIM),ZBOT(NZDIM),
     1LAYCON(NLAY),XMAX(NCOL),YMAX(NROW),NCON(NLAY),QSS(NCOL,NROW,NLAY),
     2POR(NCOL,NROW,NLPOR),QX(NCP1,NROW,NLAY),QY(NCOL,NRP1,NLAY),
     3QZ(NCOL,NROW,NLP1),DELC(NROW),DELR(NCOL),HEAD(NCOL,NROW,NLAY),
     4QSTO(NCOL,NROW,NLAY),INILOC(NPART)
C
C  THE VARIABLE "ISS" IS A FLAG INDICATING STEADY STATE OR TRANSIENT
C  CONDITIONS. IN THIS PROGRAM, ISS IS SET EQUAL TO 1 FOR STEADY STATE.
C
      KOLD=KP
      NSEGS= NCOL*NROW*NLAY
      IDSCH=1
      IBND=IBOUND(JP,IP,KP)
      HED=HEAD(JP,IP,KP)
C
C  RETURN IF PARTICLE IS IN INACTIVE CELL
C
      IF(IBND.EQ.0.OR.HED.EQ.HDRY) THEN
      IDSCH= -1
      RETURN
      END IF
C
C  RETURN IF PARTICLE IS IN THE AUTOMATIC DISCHARGE ZONE
C
      CALL ISIZS(IBOUND(JP,IP,KP),IZSTOP,IZSYES)
      IF(IZSYES.EQ.1) THEN
        IDSCH= -2
        RETURN
      END IF
C
C  INITIALIZE VARIABLES AND COMPUTE INITIAL VALUE OF
C  NORMALIZED Z COORDINATE
C
      IF(IGRID.EQ.1) THEN
      NK=KP
      NKP=NK+1
      ELSE
      NK = (KP-1)*NCOL*NROW + (IP-1)*NCOL + JP
      NKP= NK + (NCOL*NROW)
      END IF
      IF(KP.EQ.NLAY.AND.ZLOC.LT.0.0) THEN
      WRITE(I7,*)' RUN STOPPED. INCONSISTENT CONFINING BED LOCATION (1)'
      WRITE(I7,*) IPART,ZLOC,NSTEP
		call stopfile  ! emrl
      STOP
      END IF
      ZMX=ZTOP(NK)
      HED=HEAD(JP,IP,KP)
      IF(LAYCON(KP).EQ.1) ZMX=HED
      IF(LAYCON(KP).GT.1.AND.HED.LT.ZMX) ZMX=HED
      IF(ZLOC.GE.0.0) THEN
      ZP=ZLOC*ZMX + (1.0-ZLOC)*ZBOT(NK)
      ELSE
      ZL= 1.0+ZLOC
      ZP=ZL*ZBOT(NK) + (1.0-ZL)*ZTOP(NKP)
      END IF
C
C  WRITE STARTING COORDINATES
C
      IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XP,YP,ZLOC,ZP,TIME,
     1                          JP,IP,KP,NSTEP,NROW,NCOL)
C
C  CHECK TO SEE IF PARTICLE IS LOCATED IN CONFINING BED.
C  IF IT IS, MOVE PARTICLE TO NEXT ACTIVE LAYER.
C
      IF(ZLOC.GT.0.0) GO TO 50
      IF(ZLOC.GE.0.0.AND.QZ(JP,IP,KP+1).GE.0.0) GO TO 50
      IF(KP.EQ.NLAY) GO TO 30
      IF(ISS.EQ.0 .AND. NCON(KP).EQ.0) GO TO 50
      IF(NCON(KP).EQ.0) GO TO 30
      LCON=NCON(KP)
      KPOR=NLAY+LCON
      V= QZ(JP,IP,KP+1)/POR(JP,IP,KPOR)/DELC(IP)/DELR(JP)
      IF(V.EQ.0.0) THEN
      IDSCH=0
      RETURN
      END IF
      ZCB= ZTOP(NKP)
      ZMN= ZBOT(NK)
      IF(V.GT.0.0) DT= -ZLOC*(ZMN-ZCB)/V
      IF(V.LT.0.0) DT= -(1.0+ZLOC)*(ZMN-ZCB)/V
      TIM=TIME+DT
C
C  CHECK TO SEE IF MAXIMUM TIME HAS BEEN REACHED OR EXCEEDED.
C  IF SO, SET TIME EQUAL TO TMAX, CALUCULATE PARTICLE LOCATION, &
C  EXIT SUBROUTINE.
C
      IF(TMAX.GT.TIM) GO TO 10
      DTT=TMAX-TIME
      ZLC= ZLOC + V*DTT/(ZMN-ZCB)
      ZP= -ZLC*ZCB + (1.0+ZLC)*ZMN
      TTM= TMAX
      IF(ICASE.GT.1) TTM= -TMAX
      IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XP,YP,ZLC,ZP,TTM,
     1                          JP,IP,KP,NSTEP,NROW,NCOL)
      TIME=TMAX
      ZLOC=ZLC
      RETURN
10    CONTINUE
C
      TIME=TIME+DT
      ZP=ZMN
      ZLOC=0.0
      IF(V.GE.0.0) GO TO 20
      KP=KP+1
      ZP=ZCB
      ZLOC=1.0
20    IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XP,YP,ZLOC,ZP,TIME,
     1                          JP,IP,KP,NSTEP,NROW,NCOL)
      GO TO 40
30    WRITE(I7,*)' RUN STOPPED. INCONSISTENT CONFINING BED LOCATION.(2)'
      WRITE(I7,*) 'PARTICLE,  LOCAL Z,  TIME STEP,  ROW  COLUMN  LAYER'
      WRITE(I7,*) IPART,ZLOC,NSTEP,IP,JP,KP
		call stopfile  ! emrl
      STOP
40    CONTINUE
C
C  THE PARTICLE IS NOW IN AN ACTIVE MODEL LAYER.  CHECK TO SEE IF IT IS
C  IN AN INACTIVE CELL WITHIN THE MODEL LAYER.  IF SO, RESET
C  KP AND ZLOC TO THE VALUES AT THE EXIT POINT OF PREVIOUS ACTIVE CELL.
C
      IBND=IBOUND(JP,IP,KP)
      HED=HEAD(JP,IP,KP)
      IF(IBND.EQ.0.OR.HED.EQ.HDRY) THEN
      IDSCH=0
      KP=KOLD
      ZLOC= -1.0
      RETURN
      END IF
50    CONTINUE
C
C       START TRACKING LOOP
C
      DO 150 N=1,NSEGS
C
C  STORE CURRENT CELL LOCATION AND LOCAL Z COORDINATE
C
      JOLD=JP
      IOLD=IP
      KOLD=KP
C
C  ASSIGN SCALAR VARIABLES FOR CONVENIENCE
C
      XMN=0.0
      IF(JP.GT.1) XMN=XMAX(JP-1)
      XMX=XMAX(JP)
      YMN=0.0
      IF(IP.LT.NROW) YMN=YMAX(IP+1)
      YMX=YMAX(IP)
      IF(IGRID.EQ.1) THEN
      NK=KP
      ELSE
      NK= (KP-1)*NCOL*NROW + (IP-1)*NCOL + JP
      END IF
      ZMX=ZTOP(NK)
      HED=HEAD(JP,IP,KP)
      IF(LAYCON(KP).EQ.1) ZMX=HED
      IF(LAYCON(KP).GT.1.AND.HED.LT.ZMX) ZMX=HED
      ZMN= ZBOT(NK)
      IF(ZMX.EQ.ZMN) ZMX=ZMX+0.01
      ZP= ZLOC*ZMX + (1.0-ZLOC)*ZMN
      DZ= ZMX-ZMN
C
C  ASSIGN FACE VELOCITIES
C
      PHI=POR(JP,IP,KP)
      VX1= QX(JP,IP,KP)/PHI/DELC(IP)/DZ
      VX2= QX(JP+1,IP,KP)/PHI/DELC(IP)/DZ
      VY1= QY(JP,IP+1,KP)/PHI/DELR(JP)/DZ
      VY2= QY(JP,IP,KP)/PHI/DELR(JP)/DZ
      VZ1= QZ(JP,IP,KP+1)/PHI/DELC(IP)/DELR(JP)
      VZ2= QZ(JP,IP,KP)/PHI/DELC(IP)/DELR(JP)
C
C  COMPUTE CELL TRANSIT TIMES IN X, Y, AND Z DIRECTIONS
C
      CALL DTCALC (VX1,VX2,VX,DVXDX,XMN,XMX,XP,DTX,IVXFLG)
      CALL DTCALC (VY1,VY2,VY,DVYDY,YMN,YMX,YP,DTY,IVYFLG)
      CALL DTCALC (VZ1,VZ2,VZ,DVZDZ,ZMN,ZMX,ZP,DTZ,IVZFLG)
C
C  DETERMINE EXIT FACE AND SET FLAGS
C
      IX=0
      IY=0
      IZ=0
      DT=DTX
      IX=1
      IF(DTY.GE.DT) GO TO 60
      DT=DTY
      IX=0
      IY=1
      IZ=0
60    IF(DTZ.GE.DT) GO TO 70
      DT=DTZ
      IY=0
      IX=0
      IZ=1
70    CONTINUE
      IF(VX.LT.0.0.AND.IX.GT.0) IX=-IX
      IF(VY.LT.0.0.AND.IY.GT.0) IY=-IY
      IF(VZ.LT.0.0.AND.IZ.GT.0) IZ=-IZ
C
C  CHECK TO SEE IF THIS CELL IS A DISCHARGE POINT
C
C    CHECK TO SEE IF THERE IS NO OUTFLOW FROM THIS CELL.
C    IF SIMULATION IS STEADY STATE, DISCHARGE PARTICLE AND RETURN.
C
      IXYZ=0
      IF(IVXFLG.GT.1) IXYZ=IXYZ+1
      IF(IVYFLG.GT.1) IXYZ=IXYZ+1
      IF(IVZFLG.GT.1) IXYZ=IXYZ+1
 
C... IF THERE IS NO WAY OUT OF THIS CELL, START CHECKING TO SEE IF PARTICLE
C    SHOULD BE STOPPED
      IF(IXYZ.EQ.3) THEN
C
C... STEADY STATE RUNS, STOP & RETURN IF PARTICLE ENTERS CELL WITH NO
C    WAY OUT
        IF(ISS.EQ.1) THEN
          CALL ISIZS(IBOUND(JP,IP,KP),IZSTOP,IZSYES)
          IF(IZSYES.EQ.1) THEN
            IDSCH= -2
          ELSE
            IDSCH=0
          END IF
          RETURN
C
C... TRANSIENT RUNS, STOP & RETURN IF PARTICLE ENTERS CELL WITH NO WAY
C    OUT AND THERE IS A NET INTERNAL SINK (FORWARD RUN) OR NET INTERNAL
C    SOURCE (BACKWARD RUN)
        ELSE IF(ISS.EQ.0) THEN
C... IF TRACKING IS BACKWARD AND CELL HAS A NET INTERNAL SOURCE, STOP PARTICLE.
          IF(IREV.EQ.1 .AND. QSS(JP,IP,KP).GT.0.0) THEN
            CALL ISIZS(IBOUND(JP,IP,KP),IZSTOP,IZSYES)
            IF(IZSYES.EQ.1) THEN
              IDSCH= -2
            ELSE
              IDSCH=0
            END IF
            RETURN
C... IF TRACKING IS FORWARD AND CELL HAS A NET INTERNAL SINK, STOP PARTICLE
          ELSE IF(IREV.EQ.0 .AND. QSS(JP,IP,KP).LT.0.0) THEN
            CALL ISIZS(IBOUND(JP,IP,KP),IZSTOP,IZSYES)
            IF(IZSYES.EQ.1) THEN
              IDSCH= -2
            ELSE
              IDSCH=0
            END IF
            RETURN
          END IF
        END IF
      END IF
 
 
C
C...  STOP & RETURN IF PARTICLE IS IN THE AUTOMATIC DISCHARGE ZONE
      CALL ISIZS(IBOUND(JP,IP,KP),IZSTOP,IZSYES)
      IF(IZSYES.EQ.1) THEN
        IDSCH= -2
        RETURN
      END IF
 
C... KEEP TRACK OF WHETHER THE PARTICLE IS STILL IN ITS INITIAL CELL
      IF(INILOC(IPART).GT.0) THEN
      NODNUM= (KP-1)*NCOL*NROW + (IP-1)*NCOL + JP
      IF(NODNUM.NE.INILOC(IPART)) INILOC(IPART)= -INILOC(IPART)
      END IF
C
      IF(IBOUND(JP,IP,KP).GT.-1000.AND.IBOUND(JP,IP,KP).LT.1000)
     1  GO TO 80
      IBD=0
      IF(IBOUND(JP,IP,KP).LE.-1000.OR.IBOUND(JP,IP,KP).GE.1000) IBD=1
 
C... START CHECKING FOR WEAK SINK CONDITIONS:
C... IF THERE IS A SINK IN THIS CELL AND ONE OF THE TWO "STOP AT WEAK SINKS"
C    OPTIONS IS ON, THEN START CHECKING TO SEE IF THE PARTICLE SHOULD STOP.
      IF(IBD.EQ.1.AND.ISNK.GE.1) THEN
 
C... IF THIS IS THE FIRST PASS THROUGH AT TIME=0 AND THE RUN IS BACKTRACKING,
C    DON'T LET THE PARTICLE STOP IN THIS CELL
      IF(TIME.LE.0.0.AND.IREV.EQ.1) GO TO 80
 
C... IF PARTICLE IS STILL IN ITS INITIAL CELL AND RUN IS BACKTRACKING,
C    DON'T LET THE PARTICLE STOP IN THIS CELL
      IF(INILOC(IPART).GT.0 .AND. IREV.EQ.1) GO TO 80
 
      IF (ISNK.EQ.1) THEN
        CALL ISIZS(IBOUND(JP,IP,KP),IZSTOP,IZSYES)
        IF(IZSYES.EQ.1) THEN
          IDSCH= -2
        ELSE
          IDSCH=0
        END IF
        RETURN
      END IF
      END IF
C
C  DO A SECOND CHECK TO SEE IF PARTICLE CANNOT GET OUT OF CELL.
C  IF IT CANNOT, DISCHARGE IT.
C
      IF(TIME.LE.0.0.AND.IREV.EQ.1) GO TO 80
      IF(INILOC(IPART).GT.0 .AND. IREV.EQ.1) GO TO 80
      DXDY= DELR(JP)*DELC(IP)
      DXDZ= DELR(JP)*(ZMX-ZMN)
      DYDZ= DELC(IP)*(ZMX-ZMN)
      VSIGN=1.0
      IF(IREV.EQ.1) VSIGN= -1.0
      VLX1=VSIGN*VX1
      VLX2=VSIGN*VX2
      VLY1=VSIGN*VY1
      VLY2=VSIGN*VY2
      VLZ1=VSIGN*VZ1
      VLZ2=VSIGN*VZ2
      IF(QSS(JP,IP,KP).GE.0.0) GO TO 80
      QI=QSTO(JP,IP,KP)
      IF(VLX1.GT.0.0) QI=QI+(VLX1*DYDZ*POR(JP,IP,KP))
      IF(VLX2.LT.0.0) QI=QI-(VLX2*DYDZ*POR(JP,IP,KP))
      IF(VLY1.GT.0.0) QI=QI+(VLY1*DXDZ*POR(JP,IP,KP))
      IF(VLY2.LT.0.0) QI=QI-(VLY2*DXDZ*POR(JP,IP,KP))
      IF(VLZ1.GT.0.0) QI=QI+(VLZ1*DXDY*POR(JP,IP,KP))
      IF(VLZ2.LT.0.0) QI=QI-(VLZ2*DXDY*POR(JP,IP,KP))
      IF(QI.LE.0.0) GO TO 80
      F= -QSS(JP,IP,KP)/QI
      IF(F.GT.FRAC) THEN
        CALL ISIZS(IBOUND(JP,IP,KP),IZSTOP,IZSYES)
        IF(IZSYES.EQ.1) THEN
          IDSCH= -2
        ELSE
          IDSCH=0
        END IF
      RETURN
      END IF
C
80    CONTINUE
C
C  CHECK TO SEE IF TIME IS GREATER THAN OR EQUAL TO TMAX.
C  IF SO, SET TIME EQUAL TO TMAX, CALCULATE LOCATION & RETURN
C
      TIM=TIME+DT
      IF(TMAX.GT.TIM) GO TO 90
      DTT= TMAX-TIME
      CALL NEWXYZ (VX,DVXDX,VX1,VX2,DTT,XP,XMN,XMX,XPP,IVXFLG,0)
      CALL NEWXYZ (VY,DVYDY,VY1,VY2,DTT,YP,YMN,YMX,YPP,IVYFLG,0)
      CALL NEWXYZ (VZ,DVZDZ,VZ1,VZ2,DTT,ZP,ZMN,ZMX,ZPP,IVZFLG,0)
C
      ZLOC= (ZPP-ZMN)/(ZMX-ZMN)
C
C  WRITE COORINATES OF INTERMEDIATE POINTS
C
      TTM= TMAX
      IF(ICASE.GT.1) TTM= -TMAX
      IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XPP,YPP,ZLOC,ZPP,
     1                          TTM,JP,IP,KP,NSTEP,NROW,NCOL)
      XP=XPP
      YP=YPP
      ZP=ZPP
      TIME=TMAX
      RETURN
90    CONTINUE
C
C  COMPUTE PARTICLE COORDINATES AT THE EXIT POINT OF CELL (JP,IP,KP)
C
      CALL NEWXYZ (VX,DVXDX,VX1,VX2,DT,XP,XMN,XMX,XNEW,IVXFLG,IX)
      CALL NEWXYZ (VY,DVYDY,VY1,VY2,DT,YP,YMN,YMX,YNEW,IVYFLG,IY)
      CALL NEWXYZ (VZ,DVZDZ,VZ1,VZ2,DT,ZP,ZMN,ZMX,ZNEW,IVZFLG,IZ)
      XP=XNEW
      YP=YNEW
      ZP=ZNEW
      TIME=TIME+DT
C
      ZLOC= (ZP-ZMN)/(ZMX-ZMN)
      IF(ZLOC.GT.1.0) ZLOC=1.0
      IF(ZLOC.LT.0.0) ZLOC=0.0
C
      IF(IX.LT.0) JP=JP-1
      IF(IX.GT.0) JP=JP+1
      IF(IY.LT.0) IP=IP+1
      IF(IY.GT.0) IP=IP-1
      ICB=0
      IF(IZ.LT.0.AND.NCON(KP).NE.0) ICB=1
      IF(IZ.LT.0.AND.NCON(KP).EQ.0) KP=KP+1
      IF(IZ.GT.0) KP=KP-1
C
C  DISCHARGE PARTICLE IF IT REACHES BOUNDARY OF THE ACTIVE GRID
C
      IEXIT=0
      IF(JP.LT.1.OR.JP.GT.NCOL) IEXIT=1
      IF(IP.LT.1.OR.IP.GT.NROW) IEXIT=1
      IF(KP.LT.1.OR.KP.GT.NLAY) IEXIT=1
      IF(IEXIT.EQ.1) GO TO 100
      IBND=IBOUND(JP,IP,KP)
      HED=HEAD(JP,IP,KP)
      IF(IBND.EQ.0.OR.HED.EQ.HDRY) IEXIT=1
100   CONTINUE
      IF(IEXIT.EQ.1) THEN
      IDSCH=0
      JP=JOLD
      IP=IOLD
      KP=KOLD
      IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XP,YP,ZLOC,ZP,TIME,
     1                          JP,IP,KP,NSTEP,NROW,NCOL)
      RETURN
      END IF
C
C  SWITCH ZLOC FROM 0 TO 1 OR 1 TO 0 IF PARTICLE CHANGES LAYERS
C
      IF(IZ.LT.0.AND.NCON(KOLD).EQ.0) ZLOC=1.0
      IF(IZ.GT.0.AND.NCON(KP).EQ.0) ZLOC=0.0
      IF(IZ.GT.0.AND.NCON(KP).GT.0) THEN
      ZLOC= -1.0
      ICB=1
      END IF
C
C  WRITE COORDINATES OF EXIT POINT
C
      IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XP,YP,ZLOC,ZP,TIME,
     1                          JP,IP,KP,NSTEP,NROW,NCOL)
C
C  CHECK TO SEE IF PARTICLE HAS ENTERED CONFINING BED.
C  IF SO, MOVE PARTICLE THROUGH CONFINING BED TO AN ACTIVE MODEL LAYER.
C  OTHERWISE, GO TO END OF CELL LOOP.
C
      IF(ICB.EQ.0) GO TO 150
      KOLD=KP
      LCON=NCON(KP)
      KPOR=NLAY+LCON
      IF(LCON.EQ.0) GO TO 130
      IF(KP.EQ.NLAY) GO TO 130
      V= QZ(JP,IP,KP+1)/POR(JP,IP,KPOR)/DELC(IP)/DELR(JP)
      IF(V.EQ.0.0) THEN
      IDSCH=0
      RETURN
      END IF
      IF(IGRID.EQ.1) THEN
      NK= KP
      NKP=NK+1
      ELSE
      NK= (KP-1)*NCOL*NROW + (IP-1)*NCOL + JP
      NKP=NK + (NCOL*NROW)
      END IF
      ZCB= ZTOP(NKP)
      ZMN= ZBOT(NK)
      IF(V.GT.0.0) DT= -ZLOC*(ZMN-ZCB)/V
      IF(V.LT.0.0) DT= -(1.0+ZLOC)*(ZMN-ZCB)/V
C
C  CHECK TO SEE IF TIME IS GREATER THAN OR EQUAL TO TMAX.
C  IF SO, SET TIME EQUAL TO TMAX, CALCULATE LOCATION, &
C  EXIT SUBROUTINE.
C
      TIM=TIME+DT
      IF(TMAX.GT.TIM) GO TO 110
      DTT=TMAX-TIME
      ZLC= ZLOC + V*DTT/(ZMN-ZCB)
      ZP= -ZLC*ZCB + (1.0+ZLC)*ZMN
      TTM= TMAX
      IF(ICASE.GT.1) TTM= -TMAX
      IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XP,YP,ZLC,ZP,TTM,
     1                          JP,IP,KP,NSTEP,NROW,NCOL)
      TIME=TMAX
      ZLOC=ZLC
      RETURN
110   CONTINUE
C
      TIME=TIME+DT
      ZP=ZMN
      ZLOC=0.0
      IF(V.GE.0.0) GO TO 120
      KP=KP+1
      ZP=ZCB
      ZLOC=1.0
120   IF(MODE.EQ.1) CALL WRITPL(I2,ICMPCT,IPART,XP,YP,ZLOC,ZP,TIME,
     1                          JP,IP,KP,NSTEP,NROW,NCOL)
      GO TO 140
130   WRITE(I7,*)' RUN STOPPED. INCONSISTENT CONFINING BED LOCATION.(3)'
      write(i7,*) ipart,zloc,nstep
		call stopfile  ! emrl
      STOP
140   CONTINUE
C
C  THE PARTICLE IS NOW IN AN ACTIVE MODEL LAYER.  CHECK TO SEE IF IT IS
C  IN AN INACTIVE CELL WITHIN THE MODEL LAYER.  IF SO, RESET
C  KP AND ZLOC TO THE VALUES AT THE EXIT POINT OF PREVIOUS ACTIVE CELL.
C
      IBND=IBOUND(JP,IP,KP)
      HED=HEAD(JP,IP,KP)
      IF(IBND.EQ.0.OR.HED.EQ.HDRY) THEN
      IDSCH=0
      KP=KOLD
      ZLOC= -1.0
      RETURN
      END IF
150   CONTINUE
C
      RETURN
160   FORMAT(I5,1X,5(E20.12,1X),3(I3,1X),I4)
      END
 
C***** SUBROUTINE *****
      SUBROUTINE NEWXYZ (V,DVDX,V1,V2,DT,XP,XMIN,XMAX,XNEW,
     1  IVFLG,ISIDE)
C
      IF(IVFLG.EQ.1) GO TO 10
      IF(IVFLG.EQ.2) GO TO 20
      XNEW= XP + V*(EXP(DVDX*DT) - 1.0E+0)/DVDX
      GO TO 30
10    XNEW= XP + V1*DT
      GO TO 30
20    XNEW= XP
      GO TO 40
30    IF(ISIDE.GT.0) XNEW=XMAX
      IF(ISIDE.LT.0) XNEW=XMIN
40    CONTINUE
      IF(XNEW.LT.XMIN) XNEW=XMIN + 1.0E-3*(XMAX-XMIN)
      IF(XNEW.GT.XMAX) XNEW=XMAX - 1.0E-3*(XMAX-XMIN)
      IF(XNEW.EQ.XMIN.AND.V1.EQ.0.0E+0) XNEW=XMIN + 1.0E-3*(XMAX-XMIN)
      IF(XNEW.EQ.XMAX.AND.V2.EQ.0.0E+0) XNEW=XMAX - 1.0E-3*(XMAX-XMIN)
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE DTCALC (V1,V2,V,DVDX,X1,X2,XP,DT,IVFLG)
C
      IVFLG=0
      DT= 1.0E+20
      V2A= ABS(V2)
      V1A= ABS(V1)
      DV= V2-V1
      DVA= ABS(DV)
C
C  CHECK FOR A UNIFORM ZERO VELOCITY IN THIS DIRECTION
C  IF SO, GO TO 30 AND RETURN WITH DT SET TO 1E+20 AND IVFLG=2
C
      V=0.0
      IF(V2A.LE.1.0E-15.AND.V1A.LE.1.0E-15) GO TO 30
C
C  CHECK FOR A CONSTANT UNIFORM NON-ZERO VELOCITY IN THIS DIRECTION.
C  VARIABLE VVV CAN ONLY BE << 0 IF V1 AND V2 ARE IN THE SAME
C  DIRECTION. IF VVV < 1E-4, V1 AND V2 ARE ESSENTIALLY EQUAL. IF
C  THIS CONDITION EXISTS, GO TO 20 AND SET V=V1 AND COMPUTE TRAVEL
C  TIME USING CONSTANT VELOCITY. SET IVFLG=1.
C
      VV=V1A
      IF(V2A.GT.VV) VV=V2A
      VVV= DVA/VV
      IF(VVV.LE.1.0E-4) GO TO 20
C
C  COMPUTE VELOCITY CORRESPONDING TO PARTICLE'S POSTION
C
      DX= X2-X1
      DVDX= DV/DX
      XF= (XP-X1)/DX
      V= (1.0E+0-XF)*V1 + XF*V2
C
C  IF FLOW IS INTO CELL FROM BOTH SIDES THERE IS NO OUTFLOW, SO
C  GO TO 40 AND SET IVFLG=3 AND RETURN
C
      IF(V1.GE.0.0E+0.AND.V2.LE.0.0E+0) GO TO 40
C
C  IF FLOW IS OUT OF CELL ON BOTH SIDES, FIND LOCATION OF DIVIDE,
C  COMPUTE TRAVEL TIME TO EXIT FACE AND RETURN WITH IVFLG=0
C
      IF (V1.LE.0.0.AND.V2.GE.0.0) THEN
      IF (ABS(V).LE.0.0) THEN
      V= 1.0E-20
      IF (V2.LE.0.0) V= -V
      END IF
      END IF
      VR1= V1/V
      VR2= V2/V
      VR=VR1
      IF(VR.LE.0.0E+0) VR=VR2
      V1V2= V1*V2
      IF(V1V2.LE.0.0) GO TO 10
      IF(V.GT.0.0) VR=VR2
      IF(V.LT.0.0) VR=VR1
10    CONTINUE
      DT= ALOG(VR)/DVDX
      GO TO 50
20    CONTINUE
      IVFLG=1
      ZRO= 1.0E-15
      ZROM= -ZRO
      V=V1
      IF(V1.GT.ZRO) DT= (X2-XP)/V1
      IF(V1.LT.ZROM) DT= (X1-XP)/V1
      GO TO 50
30    IVFLG=2
      GO TO 50
40    IVFLG=3
50    CONTINUE
      RETURN
      END
 
C***** SUBROUTINE *****
      SUBROUTINE ISIZS(IBOUND,IZSTOP,IZSYES)
      IZSYES=0
      IF(IZSTOP.GT.1) THEN
        IB=IBOUND
        IF(IB.LT.0) IB= -IB
        IF(IB.GE.1000) IB=IB/1000
        IF(IB.EQ.IZSTOP) IZSYES=1
        END IF
      RETURN
      END
